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Abstract 

A hybrid-stress finite element method is developed for accurate stress and 
vibration analysis of problems in linear anisotropic elasticity. 

A modified form of the Hellinger-Reissner principle is formulated for 
dynamic analysis and an algorithm for the determination of the anisotropic elastic and 
compliance constants from experimental data is developed. These schemes have been 
implemented in a finite element program for static and dynamic analysis of linear 
anisotropic two-dimensional elasticity problems. 

Specific numerical examples are considered to verify the accuracy of the 
hybrid-stress approach and compare it with that of the standard displacement method, 
especially for highly anisotropic materials. 

It is that the hybrid-stress approach gives significantly better results than 
the displacement method. Preliminary work on extensions of this method to 
three-dimensional elasticity is discussed, and the stress shape functions necessary for 
this extension are included. 
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INTRODUCTION 


Objectives 

Stress-hybrid finite elements are developed for the study of equilibrium 
and vibration problems in two- and three-dimensional linear anisotropic elasticity. 

It has been observed that standard displacement-type finite element 
methods may produce very poor approximations of stresses, displacements, natural 
frequencies and mode shapes for strongly anisotropic elastic bodies, such as 
structures composed of single metallic crystals. The possibility of resolving these 
deficiencies by using elements based on assumed stresses is explored in this study. 

Further motivation of this work is a result of problems faced in the 
analysis of single crystal turbine blades subjected to very large centrifugal and fluid 
forces as in the fuel pump for the space shuttle main engine, where the high degree of 
anisotropy in the crystal leads to very poor results using standard displacement-type 
finite elements. 

The static problem is solved using a form of the Hellinger-Reissner 
energy functional where the displacements and stresses may be independently 
interpolated. To accommodate anisotropic materials, an algorithm for computing the 
elasticity and compliance constants from experimental data is developed. For the 
analysis of vibration problems, a special Hellinger-Reissner-type variational 
principle, valid for dynamic problems in linear elasticity, is also developed. The 
eigen pairs are extracted using a standard subroutine package. 

Some numerical examples in two-dimensional linear elasticity are chosen 
to verify the accuracy of the assumed- stress approach and compare the results 
obtained therein with those obtained using the standard displacement-type finite 
element method. Specifically, end loaded cantilever beams are chosen to simulate in a 
crude way turbine blades without unnecessary complexity in geometry or loading. 
Various degrees of anisotropy are assumed for the materials and the results obtained 
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are compared with analytic solutions whenever the latter are available. 

It is observed that there is a significant difference in the results obtained 
using conventional and assumed-stress finite elements, especially as the degree of 
anisotropy increases and the orientation of the local axes changes with respect to the 
global axes. The hybrid-stress elements behave very consistently and give good 
approximations of the stresses, natural frequencies and mode shapes, independent of 
the degree of anisotropy and element orientation while displacement elements are 
found to be sensitive to changes in material properties and element orientation. 


Historical Comments 

The mathematical, and practical aspects of assumed-displacement finite 
elements have been the subject of extensive research for many years. As a result of 
the sound theoretical fundamentals, the existence of a good mathematical basis and 
the ease of element formulation has resulted in the widespread use of such 
formulations. Indeed, most general-purpose computer programs employ the classical 
displacement approach. 

However, the displacement finite element model has some shortcomings 
which are evident in constrained media problems such as those involving 
incompressible materials and plate elements requiring only C° continuity 
interpolations (instead of the more difficult C 1 continuity interpolations that are based 
on thin plate theories). In these cases, "element locking" may occur as a result of very 
stiff solutions that arise as the constraint condition is approached. Reduced integration 
or selective-reduced integration can alleviate the problem partially, but an undesirable 
consequence may be the appearance of unwanted spurious energy modes. 

An alternative approach is to use the hybrid-stress model initiated by Pian 
[1], who independently interpolated intra-element stresses and compatible boundary 
displacements, using a variational principle akin to the principle of minimum 
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complementary energy. 

The ability to interpolate the stresses independently led to the solution of 
problems in fracture mechanics [2], thick, laminated composite materials [3], and 
constrained media (nearly incompressible materials) [4]. However, the versatility of 
independently interpolating stresses also leads to serious numerical difficulties. 

A minimum number of stress parameters is required to ensure correct 
stiffness rank and care has to be taken to make sure that the stiffness is invariant 
under simple transformations of the coordinates, and to prevent the entry of spurious 
energy modes into the element stiffness. Work on a systematic approach to the 
selection of these stress parameters has been done by Spilker, Maskeri and Kania [5] 
(for complete stress polynomials) and recently by Rubinstein, Punch and Atluri [6] 
and Punch and Atluri [7,8] using group theoretical methods to minimize the number 
of stress parameters and still satisfy rank and invariance conditions, for both two- and 
three-dimensional isoparametric elements. 

In addition, Pian and Chen [12] as well as Pian, Chen and Kang [10] 
have proposed new formulations for the hybrid- stress method, where the 
Hellinger-Reissner principle is used to generate elements in which equilibrium is not 
satisfied a-priori. Stress equilibrium is introduced in these methods by means of 
Lagrange multipliers. Another formulation, following the Hu-Washizu principle, is 
also suggested in [19]. 

Outline 

Following a brief discussion of some variational principles, the variational 
formulation for dynamic problems in linear anisotropic elasticity, based on a 
Hellinger-Reissner-type principle, is developed. 

Our third section, entitled "The Discrete Problem: Hybrid Stress 
Formulation", deals with the finite element discretization of the continuum problem, 
and the calculation of the element stiffness and mass matrices. A separate section is 
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devoted to the development of an algorithm for determining the elastic and compliance 
constants from experimental data in the form of Young's moduli and Poisson's ratios 
in different directions. 

In a further section, some numerical examples are presented and the 
hybrid-stress results are compared with those obtained using the displacement model. 
Problems involving materials with varying degrees of anisotropy are also considered 
to compare the two models. 

Conclusions and suggestions for further work, especially in three 
dimensions, are discussed in our final section. 


4 



THE CONTINUUM PROBLEM: VARIATIONAL FORMULATION 


Introduction 

In this section, following a brief description of some general variational 
principles, a modified form of the Hellinger-Reissner principle valid for dynamic 
problems in linear elasticity is developed. 

The classical stress-hybrid models for finite element analysis [1] are 
derived from the principle of minimum complementary energy for which the 
functional to be varied is given by 

Kq = J - o T S odv - J T t u ds (2.1) 

V 2 S u 

where a is the stress tensor (here a vector of stress components) which satisfies the 
equilibrium equation D T a + f = 0, f are the body forces, S is the compliance 
matrix, u are the prescribed displacements on the boundary S u , and T is the 
surface traction vector, where D T is the matrix divergence operator. 

Alternatively, the equilibrium conditions may be regarded as a constraint 
and a term of the form X T (D T a + f) can be introduced into the functional with the 
Lagrangian multipliers X being identified with the displacement field. Then, 

Uq = f o T S a dv + J X T (D t a + f) dv - J T t u ds (2.2) 
v 2 V s u 

Integrating the second term by parts, we get 

*hr = J J- a T S a dv - f o T (Du) dv - f T t (u - 0) ds (2.3) 

J V 2 J V JS u 
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where D is the differential operator associated with the strain -displacement relations 
(S O = e = D u). The functional is in fact the one from which the Hellinger-Reissner 
principle is developed. 

Pian and Chen [12] have recently proposed that the Hu-Washizu principle 
may also be used in the derivation of hybrid-stress elements. The functional to be 
varied is 


%w = f I - f T Ce- o T e + o t (Du)] dv - f T t (u - u) ds (2.4) 
•'y 2 J S 

where e is the strain tensor, and C = S' 1 is the elasticity matrix. Here, the 
constitutive relation for linear elasticity and equilibrium are not pre-supposed but fall 
out as the Euler- Lagrange equations for the functional t^hav- 

The last two variational formulations, though involving more computation 
time, have an advantage over the classical formulation in that equilibrium is not 
necessarily satisfied throughout the domain but only at selected points (integration 
points) in a weighted sense. This leads to a stiffness matrix that is less sensitive to 
coordinate transformation. 


Variational Formulation for the Dynamic Problem 

The variational principle that is developed is a modified form of the 
Hellinger-Reissner functional and is given by 





(2.5) 
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for a time interval [Xj,t 2 ], where T is the prescribed traction vector on the surface 
S c , p is the mass density of the material, and u is the velocity vector. The last term 
represents the kinetic energy of the body. 

Since we are primarily interested here in the eigenvalue problem, which 
is a quasi steady-state problem, the integration with time is of little consequence. 

Upon taking the first variation of the functional with respect to o and u, 

we get 


87Cmhr = 1 [ f So 8c T dv - f Du 8a T dv - [ a T D(5u) dv 
Tj V V V 

+ f (5u) T T ds - f p 8u u T dv ] dt (2.6) 

* S o * V 

Equating the first variation to zero, we obtain the following equations: 

S a = D u on V 
O • n = T on S a 

D t a - p ii T = 0 on V ^2 7 ) 

These are recognized as the constitutive relations for the elastic material, the traction 
boundary conditions, and the linear momentum equation, without body forces, 
respectively. 

For the static case, the kinetic energy term drops out and the functional 

reduces to 



( 2 . 8 ) 
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where it is understood that u = u on S u .which is similar in form to the functional 
^HR (2.3). 


THE DISCRETE PROBLEM: 
HYBRID STRESS FORMULATION 


Introduction 

In this section, we discuss the discrete hybrid-stress approximation to the 
continuous problem defined in the previous section. A detailed account of the 
calculations of the stiffness and mass matrix are given, along with derivations of the 
stress shape functions for both two- and three-dimensional finite elements. 

Following this, an extended discussion of the algorithm for the 
determination of the elasticity and compliance constants from experimental data is 
presented. 

The section ends with a note on the algorithm used to write a general 
finite element program using either displacement or hybrid-stress elements to solve 
static and dynamic problems in two-dimensional linear elasticity. 

Discrete Formulation 

The functional to be minimized in the continuum problem is given by 

ji MH r = f i a T Sa dv - J c T (Du) dv + [ u T T ds - f \ p u T u dv 
' \ 2 J v S c v 2 

(3.1) 

If the continuum is discretized into n elements then the discrete form of 
the continuous functional is given by 

n MHR = £ { ° Ts ° dv - J o T (Du) dv 
i=l «m 2 fl m 

- J - p u T u dv + J u t T ds (3.2) 

2 ^ 3£5 o m 

where £2 m represents the volume of the elements, and dC2 a represents the 

m 

boundary over which the surface traction is prescribed. 
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As shown in the previous section, the above functional (for the static 
case) reduces to the conventional assumed-stress model when the stresses are forced 
to satisfy equilibrium (in the absence of body forces). This was noted by Pian [1] 
who also demonstrated its equivalence to the Hellinger-Reissner functional with the 
stresses satisfying equilibrium. 

Thus, the stresses in each element Q m are interpolated in terms of stress 
parameters, p , in the form 

G = P p (3.3) 

where P = P(x,y,z) contains polynomial terms such that the homogeneous 
equilibrium equations 

D T G = D T (P p) = 0 (3.4) 

are satisfied exactly and identically throughout the domain 

The displacement field u is interpolated using the standard shape 
functions that are used in isoparametric elements. Hence, 

u = N($, ti, 0 q (3.5) 

(3.5) where N(£,r),£) contains the element shape functions in terms of the local 
coordinates , , and q is the vector of the element nodal displacements, such that 
the global coordinates are expressed in terms of the master element coordinates as 


x = £ Nj (5, rj, 0 Xj 
1 1 


y-f N, a, n, o yj 

z = S Nj (£, T|, Q Zj (3.6) 

where Nj are the appropriate shape functions and (x^ y i5 Zj) are the coordinates of 
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the i** 1 node, the summation being over all the nodes in an element. 

Then, the strain-displacement relation gives rise to 

e = Du = D <N(U,C) 9) = B&q.O q = y-p B* ($,q,0 q (3.7) 

where B(£,q,£) = D(N(£,q,Q) and IJI is the Jacobian of the transformation from 
local to global coordinates. 

Substituting equations (3.5), (3.3) and (3.7) into equation (3.1), and after 
noting that dv = IJI dt, dq d£ and defining the following matrices : 


and 


We obtain 


J A A 

H = J J J P T S P IJI d£ dq dC 
J 0 0 0 

(3.8) 

G =/ J 1 f 1 P T B* d£ dq dC 
0 0 0 

(3.9) 

Q =^N T Tds 

(3.10) 

M = f f 1 f 1 N T Np IJI d^ dq d£ 
0 *0 J 0 

(3.11) 

*MHR.= 1, ( 5 P T|, P - P T Gq + q T Q • -j q T Mq ) 

(3.12) 


Calculation of Element Stiffness and Mass Matrices 

As the stresses are independent from element to element, the vector of 
stress parameters P may be eliminated at the element level by taking the first 
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variation of the functional in (3.12) and equating it to zero, thus solving for P in that 
element Taking the first variation with respect to P , we get 

5 *MHR = i, (H p 8 P T - Gq 6 p T }= 0 (3.13) 

0 .= i 

Since the P's in every element are arbitrary, 

(H p - Gq) 8p T = 0 

and P = H 1 Gq . (3.14) 

Substituting back into equation (3.12), we get 

’'MHR = r|iq T G T H'‘HH- 1 Gq-q T G T H- 1 Gq+q T Q -^q T Mq| (3.15) 

1=1 2 2 

= I {- - q T G T H 1 Gq + q T Q - - q T Mq ) (3.16) 

‘ i=l 2 2 

Equating the first variation of the above functional to zero, we obtain 
(-G T H’ 1 Gq+Q - Mq ) 8 q T = 0 (3.17) 

Thus Mq + Kq = Q (3.18) 

so that 

K = G t H' 1 G (3.19) 

is the element stiffness matrix, and M as defined in (3.11) is the element mass 
matrix. 
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Displacement and Stress Shape Functions 

The displacement functions are the same as those in the 
assumed-displacement model. For two-dimensional elements, two kinds of elements 
are chosen: the four-noded bilinear quadrilateral elements and the eight-noded 
quadratic quadrilateral elements (serendipity elements). These are shown in Figure 1. 

Since the stresses inside each element are interpolated independently of 
the displacement, a multitude of choices exists for the stress polynomials that 
constitute the P matrix. A number of factors govern the choice of these polynomials, 
however. 

In the assumed displacement model, the stiffness matrix is automatically 
free from zero-energy or kinematic deformation (rigid body) modes, and is invariant 
with respect to simple transformations of the reference/natural coordinates. 

In addition to the advantages of a hybrid-stress model, viz., accurate 
stress evaluation and a stiffness matrix that is not overly rigid , it is important that the 
inherent properties of the assumed displacement model, namely, invariance with 
respect to coordinate transformations and freedom from rigid body modes of the 
stiffness matrix are preserved when the stress polynomials are chosen, for the 
hybrid-stress model to have any real use in finite element analysis. 

To avoid the presence of kinematic deformation modes, it is necessary 
that the stiffness matrix satisfy the following condition that determines its minimum 
rank (Pian [1], and Atluri and Punch [7]). 

If the number of stress parameters per element is s , the matrix H defined 
in equation (3.8) will be positive-definite and of rank s (as the energy density 
functional is always positive-definite). The rank of the stiffness matrix is then 
determined by that of the G matrix as defined in equation (3.9). 

The order of the G matrix is s x d where d is the number of generalized 
degrees of freedom per element. If there are r rigid body modes per element, then at 
best the rank of the G matrix is given by min(s, d-r), and as a consequence the rank * 
of the stiffness matrix is given by min(s, d-r). 

Since any stiffness matrix should include all the rigid body modes of the 
element, the rank of K in equation (3.19) should be d-r. 
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For a linear quadrilateral element, the minimum number of stress 
parameters is 

s = d-r = 2 x 4 - 3 = 5 . 

After satisfying equilibrium, the 9 (3 field reduces to a seven parameter 
linear stress approximation given by 

o 4 = Pj + P 2 ti + |3 6 £ 

°T[ * P 3 + ^4 5 + P 7 ^ (3 - 23) 

<^ = p 5 + M + M 

where £ andq are the local coordinates as shown in Fig. 1. 

For a quadratic quadrilateral element with 8 nodes, the minimum number 
of stress parameters is 

s = d-r = 2x8-3 = 13. 

Choosing a complete cubic stress approximation with 30 P's, upon enforcing 
equilibrium we get the following 1 8 parameter stress field : 

= Pj + P<£ + p 2 n + P 8 ti 2 + 2p 9 $Ti + p 10 £ 2 + p 13 ^ 3 
+ 3$ lA ¥r] + 3P 15 £n 2 + P 17 ti 3 

= p 3 + + + 2Pn£n + Pl2^ 2 + PlO 1 ^ 2 + 

+ P 14 T1 3 + 3P 16 ^ 2 T1 + P 18 ^ 3 
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(3.24) 


<*fy\ = P 5 + P 7 S " M “ Pll^ 2 “ P^ 2 " 2 PlO^ - 3 Pl3^ 

- 3p 14 ^T| 2 - P 15 T| 3 - P 16 ^ 3 

To reduce the number of parameters still further, we impose the 
Beltrami-Michell compatibility conditions for isotropic materials or the stress field. 
For the more general case of plane strain, in addition to equilibrium, this requires that 
(Sokolnikoff [1 1]) 

V 2 (<j£ + Oq) = 0 (3.25) 

when there are no body forces. 

Substituting (3.24) into (3.25), eliminating any redundant P's and 
renumbering, we get a 15 p stress field given by 

a 4 = Pi + P6^ + P2^ + Pg 7 ! 2 + 2 P9^ + Pic£ 2 + Pi2^ 2 + Pi3( 3 ^V2ri 3 ) 

+ 3p 14 ^ 2 -p 15 ^ 

a T 1 = P 3 + p4^ + P7 7 ! + 2 Pn^ - Pg^ 2 + PlO^ 2 - 2 & + Pl 2 ( 3 ^ 2 - 2 ^ 3 ) 

+ P 13 q 3 + 3P 15 4 2 T1-P 1 4^ 3 

= P 5 - P?£ ” Pe^ ~ Pn^ 2 ~ P^ 2 “ 2 PlO^ - 3 Pi2^ 2r » 

-3p 13 ^ 2 -P 14 Ti 3 -p I5 ^ (3.26) 

However, as mentioned earlier, the Beltrami-Michell equations in the 
form (3.25) are valid only for isotropic materials. For the general case of anisotropy, 
the compatibility conditions as expressed in terms of the strain are given by 

V x V x E = 0 (3.27) 
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For plane strain problems, this reduces to 


+ ^y _ d^Yxy 
dy^ dxdy 


where e x , £y are the components of the strain in the x and y directions and y xy is 
the shear strain. 

For the most general case of plane anisotropy, there are six independent 
constants of elasticity, and the constitutive relation is given by: 

e x = a ll°x + a 12 a y + a 16 a xy 
£y = a 21 a x + a 22 a y + a 26 a xy 

Txy = a 61°x + a 62°y + a 66 a xy (3.29) 

where a x , Cy and a X y are the stress components in the plane and 
1 

a ll~ tf 

11 XX 


- v yx ' v xy 

a i2 £ £ 

E xx E yy 


= a 


21 


__ ^xy.x _ ^XjXy 
a 16 ” *£ Q “ a 61 

E xx ' J xy 


a 22 “ 


1 


yy 


a 26“ 


'xy,y 

E yy 


T) 


y.xy 

*xy 


-a 62 


17 


(3.30) 
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given 

= Young’s modulus in the x-direction 
Eyy = Young’s modulus in the y-direction 

v„ Y = Poisson’s ration when the tensile load is along the y-direction 
v_ v = Poisson's ration when the tensile load is along the x-direction 
T|, v _ and T|_ v v = coefficients of mutual influence of the first kind that 

*y,x xy,y 

characterize stretching due to shear stresses 
Ti.. vv and ri v _ v = coefficients of mutual influence of the second kind that 

x»xy *y»Ay 

characterize shear stresses due to normal loads. 

The matrix of elastic constants is symmetric because of the existence of an 
elastic potential which necessitates that the energy be an invariant. 

Substituting equation (3.24) i.e. the self-equilibrated stress field into 
(3.29) and thence into (3.28), we get the following stress field for a completely 
anisotropic material (in two-dimensional elasticity) 

2ao 6 

<5p - Pi + P 6 £ + M + M 2 + 2 P9^ + Pl(£ 2 + Pi2<5 2 + — — Tl 3 ) 

* a ll 


+ P 13 (3pn - 


2a 12 + a 66 


l ll 


*1 3 ) + Pi 4 ( 3 ^ 2 + - 


2a 


16 

hi 


H 3 ) - ^ Pl 5 H 3 

a ll 


3i< 3i;/:+ 2 a 1? 

Oq - P 3 + M + P 7 *i - -t:M 2 + Mu 2 - “ F) + Pi,(2^T1 + 


22 


22 


2a 26 


*22 


* 2 ) 


+ Pi2( 3 ^ 


2 _ 


2a io + a, 


12 + “66 
*22 


2a 


l 3 )+P 13 (n 3 +"-^- S 3 ) P u 4 3 


11 


22 


*22 
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(3.31) 


+ P 15 (3^+— - S 3 ) 

*22 

= P5 — P7J9 — P6 1 ! — Pn^ 2 ~ P^ 2 “ ^Pjq^ 1 ! - 3 p 12 ^ 2 n 

_ 3P 13 ^ 2 - Pj4T| 3 - P 15 ^ 3 


For the 8-noded serendipity element, Punch and Atluri [7,8] have used 
group theoretical methods to minimize the number of parameters in the stress field to 
an optimum 13, by using arguments of symmetry to remove unwanted terms from the 
cubic polynomial without affecting the rank and invariance properties of the stiffness. 
This stress field with 1 3 P's is given by 

= p, + p 3 + p 8 S + p 10 T| + p2^ 2 + p 4 ri 2 - 2p 6 ^ + 2p 7 ^ri + P 12 $ti + 3p 13 ^ 2 n 

°T1 = Pi - P3 + P11S + P 9 fi - P4S 2 + p2“n 2 - - 2 Pr 5 n + P 13 1 ! 3 + 3 p J2 ^ 2 

a ^Tl = P 5 - P<£ - Ps 1 ! + P6^ 2 + Tl 2 ) + P7 (¥ - "H 2 ) - 2p2^ - 3Pi2* 2 n - 3pi3^ri 2 

(3.32) 

Comparing this field with that obtained using compatibility, we see that 
although the stresses are complete in quadratic terms, they are all incomplete cubics. 

For three-dimensional elements, viz., 8-noded and 20-noded bricks, the 
algebra involved in reducing the number of stress parameters is very tedious. 
Moreover, as shown, the compatibility conditions would change for different material 
nodes (anisotropy). Rubinstein, Punch and Atluri [6] have arrived at stress fields 
with just enough P's to ensure sufficiency for rank for the stiffness matrix. 

For 8-noded bricks, shown in Fig. 2 (a) 



The stress approximation using a self-equilibrated field is given by: 

a $ = Pi + p3 + M + MS 

°n = Pi + P 4 + “ Pi7^ ~ P 3 

°C = Pl"P4 + 2 P9C + Pl6^ 

= P2C + P5C “ p7^1 - p8^ + Pjo + Pj 3 ^l + Pl 4 ^ 
a 4C = p 2 ri + p 6 T| - p«£ + p n - p 13 C + p 15 S 

a qC = p2^ ~ Ps^ “ Ps^ - PsC - p9^ + Pl2 " P 14 C " P^ (3.33) 

For 20-node bricks, shown in Fig. 2 (b) 
s min = d-r = 20 x 3 - 8 = 54 

Following Rubinstein, Atluri and Punch [6], the stress field that satisfies the 
equilibrium conditions is given by 

= Pi+ P 3 + m + w + 2 Pio^ + p| 6 i + m 3 +i 5 23( i i 2 + ; 2 > + iw-Ti 2 - ? 2 > 

+ 2 P 25 C 1 2 + C 2 ) - Mn 2 + ; 2 ) + 2p 27 (n 2 - C 2 > + Pas (C 2 - n 2 ) 

+ Pa&i + 2p J0 « - 6 P 4J «ti + P 46 nC 2 + P 4S Cn 2 

“ Pl~ P3 + P 4 + P 8 £ + P9C + 2Pj ]T1 + Pj 7 £ “ P]gC + P 7 2 t 1 2 + ?23(C 2 + £ 2 ) 

- fe 4 « 2 - t 2 ) - p2 5 « 2 + i 2 ) + 2P 26 (^ 2 + C 2 ) + p27« 2 - S 2 ) 

- 2p 29 $r, + 2pj,nc + p 42 « + «p 44 5 <ic + «p 45 ^; + p 47 « 2 + p 48 « 2 

= Pl — P 4 + P 7 1 + P 8 £ + 2 Pl2^ - Pi 6^ “ Pl7^ + P22^ 2 + P23^ 2 + 7) 2 ) 

+ P 24 ('n 2 - 5 2 ) - p 2 sft 2 + n 2 ) - Pj«ft 2 + 1 2 ) + p 2 7 (Ti 2 - K 2 ) 

+ Pas 01 2 - & - 2P 30 « - 2p 31 nc + p 41 4n - ep^nC + P 46 ^ 2 n 
+ P 4 7^n 2 
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- ? 2 ^ + PsC “ PlO 1 ! “ Pll£ + P 13 + Pi? 1 ! + P 20 ^ " P 22 ^ 1 ! + ?29^P~P^ 

+ 2 P 32 k - 2 p 33 nc + p^-n 2 ) + 2 p 35 « + 2 p 36 nC + p 40 c 2 

- PW^+tl 2 ) + p 45 ?(n 2 - 5 2 ) + P 50 OI 3 - 3riC 2 ) + P 5 ,(5 3 -3K 2 ) 

- P 5 3 ( 1 3 + JlC 2 ) + P(5 3 +HC 2 ) 

o 5C = P 2 n + P 6 n - Pi 0 t - p, 2 * + p, 4 - P 19 c +p 2 ,5 - P 22 « +p 30 (5 2 - ; 2 ) 

- 2p3^n + P 33 (C 2 - P) + 2p 34 nC + 2 p 35 iri - p 36 (S 2 + p) + 2 p„nt 
+ P 39 n 2 + P 44 n<2^ 2 + c 2 ) +P 45 n(s 2 + 2 p) + p 49 (S 3 - 3^n 2 ) 

+ p 50 (C 3 - 3n 2 C) - P 52 (^ 3 +3n 2 4) + P 53 (C 3 +3n 2 0 


- ?2^ PjS _ “Pi lC - P12I + P15 “ P 2 0^ “ P2I 1 ! “ P22 1 !^ 

+ p 3I (n 2 - C 2 ) + P 32 (n 2 - ; 2 ) + 2p 33 5n - 2 p 34 « - p 35 <n 2 + p> 
+2p 36 ^Ti + 2p„« + p 38 52 +P 44 ^(t] 2 ~C 2 ) - P^Ol 2 +2t 2 ) 

+ P^tt 3 - 3^ 2 n) P 5 ,« 3 - 3pO + P 52 (n 3 + H 2 ^ - p 54 (C 3 + 3 | 2 o 

(3.34) 


Elastic Constants from Experimental Data 

For an isotropic medium, the measurement of the elastic constants is a 
straightforward process. The measurements can be made along the specimen axes 
without regard to the crystal structure of the material, as the properties are the same in 
all directions. By contrast, there could be as many as 21 independent constants of 
elasticity in a fully anisotropic material instead of just 2 as in an isotropic material. 
Measurement of the elastic (Young's) moduli and shear moduli is a standard process, 
but direct measurement of the off-diagonal terms such as the coefficients of mutual 
influence, the coefficients of Chentsov and Poisson's ratio is difficult, to say the 
least. Thus, to determine these coefficients it is necessary to take advantage of their 
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interrelationship under rotation. 

If the measurements of the material properties are taken along an 
orthogonal coordinate system which is rotated with respect to the primary material 
axes, the compliance (or elastic) matrix with respect to the rotated coordinate system 
will be a transformation of the matrix with respect to the original coordinate system 
dependent solely on the direction cosines of the rotation of the coordinate system 
from the primary material axes to the global coordinate axes. 

In the realm of linear isothermal elasticity, there are a maximum of 21 
independent constants. Hence, any quantity that is measured in the rotated (global) 
frame of reference is a function of all 21 components of the compliance matrix in the 
unrotated frame. For a unique determination of the coefficients of the unrotated 
compliance matrix, 21 independent equations in 21 unknowns are required. These 
could be obtained by the measurement of the same quantity in 21 independent 
directions or by measuring more than one quantity in many independent directions so 
that there are at least 21 independent equations. Thus, a least square fit can be carried 
out on the data to determine the best approximation to the elastic constants in the 
compliance matrix. 

Using indicial notation in a Cartesian coordinate system, the position of 
the rotated coordinate axis with respect to the unrotated axis is given by the direction 
cosine matrix as shown in Table 3.1. 


Table 3.1 Direction Cosines 
\ x y z 


*il 

X-i2 

X.i3 

^21 

hi 

^23 

hi 

hi 

*33 


Letting x', y', z' refer to the rotated and x, y, z refer to the unrotated axes, then Xjj 
denotes the direction cosine of the angle between the xj axis and the xj axis. 
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Then, the stresses in the rotated frame of reference are given by: 


2 2 2 

a ll' = a ll^ll + °22^12 + a 33 ^13 + 2x 12^1 1^12 + 2x 13 ^11^13 + 2x 23 ^12^13 
2 2 2 

°22 = a 11^21 + °22^22 + a 33 ^23 + 2x 12^2 1^22 + 2x 13^2 1^23 + 2x 23^22^23 
2 2 2 

a 33' = a 11^31 + °22^32 + a 33 ^33 + 2x 12^31^32 + 2x 13^31^33 + 2x 23^32^33 

x 12* = CT 1 1^1 1^21 + a 22^12^22 + a 33^13^23 + x 12^11^22 + ^21^12^ 

+ + X^jX^) + T 2 3 (^12^23 + ^22^13^ 

x 13' = °1 1^1 1^31 + a 22^12^32 + °33^13^33 + x 12^il^22 + ^21^12^ 

+ ^(X^X^ + X 21 X 13 ) + x 23 (Xi 2 X 23 +X 21 X 13 ) 

x 23' = CT 1 1^21^31 + °22^22^32 + °33^23^33 + x 12^1^32 + ^31^22^ 

+ ^^(Xj jX 33 + X 3 iX 23 ) + x 23 (X 22 X 33 +X 32 X 23 ) 

Similarly, the stresses in the unrotated frame may be expressed in terms 
of the stresses in the rotated frame as 

• 2 • 2 ’ • 2 ' ’ ’ 

°11 = c ll^ll + °22^21 + °33 ^31 + 2 x 12^11^21 + 2x 13 ^21^31 + 2x 23 ^11^31 

* 2 * 2 * 2 * * ' 

°22 = C 1 1^-12 + °22^22 + a 33 ^32 + 2x 12^ 12^22 + 2x l 3^12^32 + 2x 23^22^32 

• 2 * 2 * 2 ’ * * 

°33 = °1 1^13 + °22^23 + a 33 ^33 + 2x 12*13*23 + 2x 13^13^33 + 2x 23 X2 3 X 33 


x 12 = °1 1^1 1^21 + a 22^21^22 + a 33^31^32 + x 12^11^22 + ^21^12^ 


+ x 13 (X n X 32 + X 31 X 12 ) + x 23 (X 21 X 32 +X 31 X 22 ) 
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*13 ~ °1 1^1 1^-13 + °22^21^23 + °33^31^33 + x 12^11^23 + ^21^13^ 


I I 

+ *13(^11^-33 + ^31^13) + *23 (^21^33 + ^31^23^ 


tilt 

*23 = °1 1^12^-13 + a 22^22^23 + a 33^32^33 + x 12^12^23 + ^22^1 3^ 

I t 

+ *13(^12^33 ^32^13) *23 (^22^33 + ^32^23^ (3.36) 


Now, in linear elasticity there exists an elastic potential which is a 
quadratic form in the stress components given by 

^ ~ 2 ^jkl ®ij ^kl = 2 a ijkl ^ij °kl (3.37) 

where, as in tensor notation, all the indices range from 1 to 3. 


Instead of expressing the stresses as second order tensors and the 
elasticity constants as a fourth order tensor, we consider the stresses as vectors of 
order 6x1 and the elasticity matrix as a matraix of size 6x6. Then, 


V = 1/2 X X ay Oj Oj = 1/2 X X a'jj Oj' a- (i,j = 1 , ... , 6) 


(3.38) 


Substituting the expressions for the stresses in the unrotated frame of 
reference from equation (3.37) into (3.38) we can get equations for each of the 
compliance constants in the rotated frame, i.e. a'y, in terms of the compliance 
constants in the unrotated frame, i.e. a^, and the direction cosine of the rotation. 

Then, 


a’ij = XXamnqimqjn (m,n = 1, ... , 6) 


(3.39) 
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where , qj n are the components of the transformation matrix given by 


*11 

x 2 

a 21 

x 2 

a 31 

2*n*2i 

2*11*31 

2*21*31 

*12 

k 22 

y 2 

a 32 

2*12*22 

2*12*32 

2*22 *32 

x 2 
K 13 

X 2 

a 23 

V 2 

a 33 

2*13*23 

2A.j3X.33 

2^23^33 

*11*12 

*21*22 

*31*32 

*11*22 + *21*12 

*H*32+*3i*22 

*21*32 + *31*22 

*11*13 

*21*23 

*31*31 

*H*23+*2i*i3 

*11*33+*31*1 3 

*2i*33+*3i*23 

*12*13 

*22*23 

*32*33 

*12*23 + *22*13 

*12*33+*32*i3 

*22*33^*32*23 


(3.40) 

Hence, when the compliance constants in the unrotated frame of reference 
are known, it is possible to calculate the compliance constants in any other orthogonal 
coordinate system, given the appropriate direction cosines. 

In general, the material is defined in the coordinate system of maximum 
symmetry using a bare minimum number of required parameters. However, arbitrary 
rotations can cause the generation of off-diagonal terms in the compliance matrix even 
if they were absent in the compliance matrix as defined in the primary coordinate 
system. It is thus necessary that the analysis always have the ability to deal with a 
fully anisotropic compliance matrix. 

The experimental data given was in the form of Young's moduli and the 
shear moduli, measured in several independent directions. As the rotation angles are 
known for a given measurement, the direction cosines can be calculated and the set of 
equations given by (3.39) becomes linear in the coefficients of the unrotated 
compliance matrix. Since more than the minimum number of measurements is made, 
a least squares approximation to the compliance matrix can be made. 
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For the problem at hand, viz., the turbine blade of the fuel pump in the 
Space Shuttle Main Engine, the material is a PWA nickel alloy, each blade consisting 
of 2 or 3 crystals. These crystals exhibit cubic syngony and have only three 
independent constants, the Young's modulus (E), and the shear modulus (G) which 
are the same in all three primary directions and Poisson’s ratio (v) which is 
independent of the other two constants. The compliance matrix along the material's 
primary axes is 


a ll 

a 12 

a 12 

0 

0 

0 

a 12 

a ll 

a 12 

0 

0 

0 

a 12 

a 12 

a ll 

0 

0 

0 

0 

0 

0 

a 44 

0 

0 

0 

0 

0 

0 

a 44 

0 

0 

0 

0 

0 

0 

a 44 


where ajj = 1/E ; a 12 = - v/E ; a^ = 1/G ; but when measured in arbitrary 

directions, the cross-coupling terms no longer remain zero. Given a sufficient number 
of equations, it is then possible to calculate ajj, a 12 and a 44 using a standard least 
squares technique. 

Program Description and Algorithm 

A simple finite element code is developed to compare the performance of 
the hybrid- stress model with that of the standard displacement model. The algorithm 
detailing the flow in the program is depicted in Fig. 3. 

The program consists of four main sections: 

i) the pre-processor which reads in the data that define the finite element 
mesh and the material properties as well as the boundary conditions. The 
pre-processor also sets up the element integration point co-ordinates, the appropriate 
weights and the coefficients of the compliance matrix in the global co-ordinate 
system. 
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ii) the processor which calculates the individual element stiffness matrices 
as in equation (3.19) and assembles them to form a global stiffness matrix. It also 
calculates the individual force vectors and assembles them to form a global force 
vector. The processor then enforces the various boundary conditions, viz., 
prescribed displacement, surface traction and nodal forces, and finally solves for the 
unknown nodal displacements using a standard solver. 

iii) the post-processor which calculates the stresses in each element once 
the nodal displacements are known. Unlike in the standard displacement model, 
where the strains are calculated first and then the stresses are obtained by using the 
stress-strain law, in the hybrid model, the stresses are calculated directly using 
equation (3.14), i.e. 

P = H* 1 Gq (3.14) 

where the q's are the generalized nodal displacements in an element. Once the P's 
are known, the stresses are given as 

a = P p (3.3) 

and are calculated at the Gauss points of integration. 

iv) the eigensolver, which sets up the element mass matrices as in 
equation (3.11), assembles them to form a global mass matrix, and after weighting 
the diagonal terms that correspond to prescribed displacement boundary conditions, 
calls a generalized eigensolver subroutine to calculate the eigen-pairs for the specified 
problem. 

The Eigensolver used for these calculations was based on the 
Householder scheme. 

The main feature of the complete code is that even though it has been set 
up to solve 2-dimensional static and dynamic problems, it is possible to solve 
3-dimensional problems by making only minor modifications. 

The number of degrees of freedom per node may be easily changed by 
means of a parameter statement, and the addition of 3-dimensional shape functions. 
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for displacements and stresses as well as Gaussian integration can make this program 
completely general. 
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NUMERICAL RESULTS 


Introduction 

In this section, initially, a brief discussion of the standard 
displacement model is presented, since it is used as a yardstick of comparison for the 
hybrid-stress model. The primary differences between the two schemes of analysis 
are also highlighted. 

Following this, some examples for the statically loaded case are 
studied, with an emphasis on the displacements and stresses in a cantilever beam 
subjected to an end shear. Both hybrid-stress and displacement finite elements are 
used to analyze beams that are isotropic and anisotropic. 

Then, an eigenvalue analysis of the cantilever beams is carried out to 
determine the first few natural frequencies and mode shapes, using both models for 
isotropic and anisotropic materials. 

All the results are compared to analytical solutions whenever the latter 

are available. 

The Assumed Displacement Approach 

The variational principle that is used in the standard displacement 
approach is akin to the principle of minimum potential energy and is of the form 

7 t PE = f 2 J ^ pu T udvdt + J 2 f ^e T Cedvdt 
Xj 2 Xj v 2 

- f 2 [ o T u ds dt (4.1) 

J x 1 J dv 

when there are no body forces, where 
e = the strain tensor 
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C = the elasticity tensor = S' 1 
u = prescribed displacement on the boundary S u 
T = prescribed traction on the boundary S a 
S u U S a = dv , the total boundary. 

On taking the first variation of the functional given in (4.1) and 
equating it to zero, we get the equations of motion and the traction boundary 
condition, provided that the variation in the displacement 8u is zero on the boundary 
S u , where the displacements are prescribed, i.e. 

pu - D T a =0 on V 

u = u on S u (4.2) 

a . n = T on S c 

If the continuum is divided into n discrete elements, then the 
discrete form of the functional 7t pE is given by 



f 2 f o T u ds dt} (4.3) 

where fi m represents the volume of the element and dQ m is the part of the 
boundary of the element that has either prescribed displacement or prescribed traction. 
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As the problem being solved is not truly transient but that of 
steady-state vibration, the integration in time may be left out without altering the 
results, so that 



( 4 . 4 ) 


Interpreting the displacements within each element by means of 
shape functions that are expressed in terms of Legendre polynomials, we have 


u = N(£, T|, Q q 


( 4 . 5 ) 


where N(£, r|, 0 contain the element shape functins in terms of the local 
co-ordinates T|, £ and q is the vector of generalized nodal displacements. 

The global co-ordinates expressed in terms of the shape functions are 


X = I Nj (£, T\, 0 Xj 
i 

y = I Nj (£, q, 0 yj (4.6) 

z = I Nj (£, q, 0 zj 
i 

the summation being carried out over all the nodes in each element. 

The strain-displacement relation then leads to 

e = Du = D[N($,TU)q] 

= B£,Ti,C)q (4.7) 
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where B(5, r\, 0 = D[N& T|, ©]. 


Substituting equations (4.5) and (4.7) into (4.4) and defining the 
following matrices 



K = f 1 / / B t CB IJI d£ dT\ dC 

0 0 J o 

(4.8) 

and 

,1 A A 

M = j J f p N T N |JI d^ dll dC 
o o t 

(4.9) 


Q = J N T T ds 

J S 

(4.10) 

we get 



^SSPE = 

+iq T Kq - q T Q} 

(4.11) 

zero, we get 

Taking the first variation of the above functional and equating it to 


n 

X {Mq + Kq - Q} 5q T = 0 

(4.12) 


Since the 5q are independent from element to element, this implies 
that 

Mq + Kq = Q (4.13) 

so that K and M as defined by equations (4.8) and (4.9) respectively are the 
element stiffness and mass matrices respectively, and Q as defined in (4.10) is the 
force vector. 
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For the case of static loading, the acceleration is zero, so that 


Kq = Q (4.14) 

while for the case of free vibration, 

Mq + Kq = 0 (4.15) 

For purposes of comparison with the hybrid stress model, two types 
of displacement elements are considered, namely, the four-noded linear quadrilateral 
and the eight-noded serendipity element. A finite element program incorporating 
these assumed displacement-type elements has been developed by modifying the 
stiffness matrix calculations as well as the stress evaluation routine. 

Unlike in the hybrid-stress model, where intra-element stresses are 
calculated without numerical differentiation of the displacements, the strains (and 
hence the stresses) in the assumed-displacement model are obtained in the following 
way: 


e = D [N(£, rj, Q] q (4.16) 

and the stresses are thus 

o = Ce (4.17) 

being calculated at the integration points in each element. 

Numerical examples are considered in the following paragraphs . 

Static Analysis 

Some problems in plane elasticity are analyzed using both the 
displacement and the hybrid models, and the results are compared. 
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Since a cantilever beam clamped at one end and in plane stress is the 
closest two dimensional analogue to the three dimensional turbine blade that is to be 
evenutally analyzed, it is thoroughly analyzed using the two models for isotropic and 
anisotropic cases. 

The clamped beam subjected to an end shear load is shown in Fig. 4. 


For the case when the beam is fully isotropic, classical beam theory 
gives the following results for the tip displacement, maximum bending stress and 
maximum shear stress: 


_ Wl 3 
Uti P " 3EI 


(4.18) 


°bending = 


Wy 


max 


I 


(4.19) 


and 


_VQ™x 

T shear “ 


(4.20) 


where following usual notation is used: 


W 

1 

E 

I 


ymax 


V 

Qmax 

b 


= end shear load 
= length of the beam 
= Young’s modulus 

= moment of inertia of the cross-section about the neutral 
(centroidal) axis. 

= distance from the neutral axis to the farthest point on the 
beam 

= half the depth of the beam 
= shear force at the particular cross section 
= first moment of the area above the neutral axis with 
respect to the neutral axis 
= width of the beam 
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In the case at hand, 


W = 250 lb. ; = 10" ; h = 1" ; E = 3 x 10 7 lb./in. 2 
For a rectangular cross section, I = 1/12 bh 3 ; b =1" 


Ujjp = 0.03333" 

^bending max. = 15000 lb./in. 2 
^shear max. ~ lb./in. 2 

Linear Element Results 

Using the displacement and the hybrid-stress models, five different 
cases were run for an isotropic beam, with the number of linear four-noded elements 
varying from 10 to 80. The various meshes used in the analysis are shown in Fig. 5. 

For comparison with analytically obtained results, the normalized tip 
displacement (Uhybj-jdAiana].) and the normalized maximum bending stress 
( a hybrid/ CT anal.) 3X6 plotted against the number of elements, and are shown in Fig. 6 
and Fig. 7. 

It is seen that even for the isotropic case, the hybrid model converges 
to the analytical solution faster than the assumed displacement model, for identical 
finite element meshes. Both displacements and stresses obtained using the 
hybrid-stress model are consistently better than those obtained using the assumed 
displacement model. 

For the anisotropic case, the material model chosen is that of cubic 
syngony to simulate the single crystal turbine blade made of the nickel alloy. 

From the experimental data supplied, the following material 
properties were obtained using the materials subroutine described in the previous 
section discussing Elastic Constants from Experimental Data: 
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E, = 1.9716 x 10 7 lb./in. 2 

E 2 = 1.9716 x 10 7 lb./in. 2 [neglecting the small difference from 

Ej due to possible errors in measurement]. 

= 0.2875 

G 12 = 5.4758 x 10 6 lb./in. 2 

Using these material properties, the same problem, viz., an 
end-loaded cantilever, is solved using hybrid and displacement methods. 

Since the mesh with 40 elements gave very good results for the 
isotropic case without taking up too much computation time, it is used. In the 
absence of an analytical solution, the comparison between the two models is made on 
the basis of a rotation of the material axes through various angles relative to the global 
axes. Since the two moduli of elasticity Ej and E 2 are identical, it is to be expected 
that any rotation by pairs of angles that are complementary produce the same result as 
the physical problem remains unchanged. On rotation of the axis coinciding with Ej 
by 30" , 45* , 60° , and 90° and comparing the results, shown in Table 4.1, it is just 
as expected. While the results for rotation of the material axes by 30° and 60° using 
hybrid elements are exactly the same, there is quite a variation in the results obtained 
using displacement elements. This is also seen for the rotation by 90° as compared to 
no rotation. 

Table 4.1 

Anisotropic cantilever solution using 40 linear elements for various 
orientations of the material axes. 


Orientation of the Axes 

u tip(displ.) 

“tipOtybrid) 

°max(displ.) 

°max(hybrid) 

0° 

0.046742" 

0.049624" 

11089 psi 

11117 psi 

30° 

0.047149" 

0.057905" 

10981 psi 

10579 psi 

45° 

0.051876" 

0.061880" 

12160 psi 

11023 psi 

60° 

0.048390" 

0.057905" 

10662 psi 

10579 psi 

90° 

0.050489" 

0.049624" 

10618 psi 

11117 psi 
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It is observed that there is an error of as much as 8-percent in the tip 
displacement, when the material axes are rotated by 90* , using the displacement 
model, while there is no change when the results of the hybrid model are compared. 
A similar error is noticed on comparing the results for 30* and 60* rotation when the 
displacement model is used. 

To study the difference between the two models further, an arbitrary 
anisotropic material model is considered where, 

Ej = 3.0 x 10 7 lb./in. 2 

E 2 = 3.0 x 10 6 lb./in. 2 

v 12 = 0-3 

G 12 = 1.5 x 10 6 lb./in 2 

The same problem is solved using both 40 and 80 linear elements for 
various orientations of the material axes, and the tip displacements are tabulated in 
Table 4.2. 


Table 4.2 

Tip Displacement Convergence on Refining the Mesh 

40 Linear Elements 80 Linear Elements 


Orientation of the Axis 

“tipOybrid) 

u tip(displ.) 

“tipOiybrid) 

u tip(displ.) 

0* 

0.03312" 

0.02970" 

0.03349" 

0.03231" 

30* 

0.03991" 

0.14505" 

0.04209" 

0.23364" 

45* 

0.07179" 

0.10085" 

0.07818" 

0.12172" 

60* 

0.14492" 

0.03075" 

0.15323" 

0.03295" 

90* 

0.18661" 

0.11522" 

0.18143" 

0.14632" 


From Table 4.2, it is observed that while the tip displacement 
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obtained using the assumed displacement method changes drastically on refining the 
mesh, the tip displacement obtained using the hybrid-stress method does not change 
much when the number of elements is increased, indicating that the hybrid model 
converges faster than the standard displacement model even for completely 
anisotropic materials. 

Another observation made is that the tip displacement continuously 
increases as the angle of rotation is changed from 0' to 90* for the hybrid-stress 
model, while it fluctuates arbitrarily for the displacement model. Since E 2 is a tenth 
of Ej, it is to be expected that the displacement increase as the rotation increases, 
reaching a maximum when the orientation of the material axes is 90’ away from the 
global axes. 

The hybrid model thus gives good results even for arbitrary 
anisotropic materials with material axes not coinciding with the global axes. 

Quadratic Element Results 

To compare the actual values of the tip displacements and bending 
stresses, 8 noded quadratic elements of the assumed-displacement and 
assumed-stress type are used to solve the same problem. The number of elements is 
varied from 3 to 20 and the different finite element meshes used are shown in Fig. 8. 

The normalized tip displacement (Uhybric^anal.) and the normalized 
bending stress (tfhybrid^anal.) are plotted against the number of elements, and are 
shown in Fig. 9 and Fig. 10. 

For the isotropic case, using identical grids, it is observed that the 
hybrid model converges to the analytical solution faster than the displacement model, 
both the tip displacement and the bending stress being consistently better for the 
hybrid model. 

The anisotropic material model chosen is the same as before, viz. a 
crystal with cubic syngony, to simulate the nickel alloy turbine blade. 

The 10 element mesh is chosen, and the results for various rotations 
of the material axes are presented in Table 4.3 which follows. 
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Table 4.3 


Anisotropic cantilever solution using 10 quadratic elements for 
various orientations of the material axes. 


Orientation of the Axes 

“tijKdispl.) 

“tipChybrid) 

°max(displ.) 

°max(hybrid) 

0* 

0.050798" 

0.051166" 

11489 psi 

11532 psi 

30* 

0.051840" 

0.060461" 

11857 psi 

12272 psi 

45* 

0.060356" 

0.063769" 

12014 psi 

12566 psi 

60* 

0.054663" 

0.060464" 

12099 psi 

12273 psi 

90* 

0.058364" 

0.051166" 

121 12 psi 

11532 psi 


It is seen that for rotations of the material axes by 30* and 60* 
respectively, the hybrid model gives the same tip displacement and maximum bending 
stress as expected since the moduli of elasticity Ej and E 2 are identical. The 
results from the displacement model however vary by about 6 percent for the tip 
displacement for the same rotations of the axes, and comparing the results between no 
rotation and a 90* rotation, the difference is observed to be as much as 15 percent. 

In Table 4.4 the values of the tip displacement for the cantilever with 
the cubic syngony material model using both 40 linear elements and 10 quadratic 
elements are shown. It is seen that while the results of the hybrid elements for both 
elements are very close, this is not true of the displacement elements. 
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Table 4.4 


Comparison of tip displacements for linear and quadratic hybrid 
and displacement finite elements 


u tip(displ.) u tip(hybrid) 


Orientation of Axis 

linear elements 

quadratic elements 

Linear elements 

quadratic elements 

O' 

0.046742" 

0.050798" 

0.049624" 

0.051166" 

o 

« 

0.047149" 

0.051840" 

0.057905" 

0.060461" 

45* 

0.051876" 

0.060356" 

0.061880" 

0.063769" 

60' 

0.048390" 

0.054663" 

0.057905" 

0.060464" 

90' 

0.050489" 

0.058364" 

0.049624" 

0.051166" 


This difference in the results of the displacement elements as compared 
to the hybrid elements is shown in Fig. 11. 

Next a cantilever beam with a highly anisotropic material is considered 
with the following material properties: 

Ej = 3 x 10 7 lb./in. 2 

E 2 = 3 x 10 5 lb./in. 2 

V 12 = 0.3 

G 12 = 1.5 x 10 7 lb./in. 2 

Two meshes are used, one with 10 elements and another with 20 
elements. The tip displacements using hybrid and displacement elements are shown 
in Table 4.5. 
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Table 4.5 


Convergence of tip displacement on refining the mesh 


10 quadratic elements 20 quadratic elements 


Orientation of Axis 

“tipOiybrid) 

u tip(displ.) 

“tipOiybrid) 

u tip(displ.) 

0* 

0.033597 

0.033487 

0.033694 

0.033593 

30' 

0.21054 

2.8188 

0.20560 

2.9842 

45' 

0.77374 

1.3933 

0.75420 

1.4500 

60' 

1.7581 

0.055362 

1.7513 

0.05893 

90' 

3.3327 

1.7271 

3.3431 

1.7744 


The displacement method gives results which are very different from 
those of the hybrid method. As with linear elements, since E 2 is one hundreth of 
Ej , the tip displacement should increase with rotation, reaching a maximum when 
O = 90' . This behavior, however, is exhibited only by the hybrid model, 
suggesting that it is far more stable under rotations for highly anisotropic materials. 


Vibration Analysis 

A cantilever beam is next analyzed for its first few natural frequencies 
(eigenvalues) and mode shapes using the assumed-displacement and the 
hybrid-stress method. 

As already discussed, a consistent mass is generated and the generalized 
eigenvalue problem 


Mq + Kq = 0 (4.21) 

is solved for its eigenpairs which arc the natural frequencies and mode shapes of the 
physical system. 
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Since the size of the matrices is not very large, a solver from IMSL that 
determines all the eigenvalues and eigenvectors is used instead of the sub-space 
iteration scheme suggested by Bathe [13]. 

For Bemoulli-Euler beams made of isotropic materials, neglecting the 
effect of shear deformation and rotatory inertia as we will consider only the first two 
modes where the correction introduced as a result of these effects is small, the 
equation of motion for transverse vibration is 


_d^ 

dx 2 


dV d^ 

(EI_) + pA — 

dx 2 dt 2 


0 


(4.22) 


where v = v(x,t) is the transverse displacement, 

A = area of cross section of the beam 
x = axial distance from the point of support 
p = mass density of the material 
I = centroidal moment of inertia of the cross-section 


The boundary conditions for the cantilever are: 


v = 0 and 


£-° 


d v 

— = 0 and 
dx^ 



At the fixed end; 


At the free end; 


(4.23) 


Substituting the boundary conditions into general solution, we get three 
homogeneous linear algebraic equations which yield a non-trivial solution only if the 
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determinant of the coefficients vanishes, i.e. 


1 + Cos XL CoshXL +1=0 (4.24) 

which is the characteristic equation whose roots are the eigenvalues X,. times length 
L. A numerical solution exists for the above equation, determined by Craig and 
Chang [14]. The first few values are 

XjL = 1.8751 
X^L = 4.6941 

and the natural frequencies for the cantilever are given by 

(*t L ) 2 /EIx 1/2 
L 2 V P A ' 

3.516 / El 
L 2 V P A 

22.03 (EI) 1/2 

L 2 ^ P A (4.25) 

Substituting the numerical values for the given problem, we get 

©! = 17.58 Hz. , to 2 = 110.15 Hz. 

The mode shapes are given by 

V r (x) = Cosh (XjX) - Cos (XjX) - k, [Sinh (X^) - Sin (X^)] (4.26) 

where 


1/2 

) 


so that 


and 


cor 


co, 


©2 
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(4.27) 


Cosh (XjL) + Cos (AjL) 

^ Sinh (XjL) + Sin (\L) 
as in Craig [20]. 

The first two mode shapes for a cantilever in free vibration are shown in 

Fig. 12. 

Linear Element Results 

The meshes that were used for the static problem are used here, with the 
number of elements varying from 10 to 80. 

The normalized natural frequencies (coj anaiVcOj f e and anal f.e.) 
are plotted against the number of elements, and are shown in Fig. 13 . 

For the isotropic case, the hybrid model converges to the analytical 
solution faster than the assumed displacement model. The mode shapes however do 
not seem to vary much, as seen in Fig. 14 (for the mesh with 80 elements). 

The material model chosen for the anisotropic case is that of cubic 
syngony with the same properties as in the static case, and the 40 element mesh. The 
first two natural frequencies, for various angles of rotation of the material axes, using 
both finite element approximations are tabulated in Table 4.6. 

Table 4.6 


Natural frequencies for an anisotropic cantilever beam 
using 40 linear elements for different material axes' orientations. 


Orientation of the Axes 

ti)j (hybrid) 

&)] (displ.) 

002 (hybrid) 

^(displ.) 

0* 

14.351 Hz 

14.785 Hz 

85.950 Hz 

88.341 Hz 

30* 

13.307 Hz 

14.726 Hz 

80.673 Hz 

88.0945 Hz 

45* 

12.873 Hz 

14.058 Hz 

78.379 Hz 

85.011 Hz 

60* 

13.307 Hz 

14.547 Hz 

80.673 Hz 

87.318 Hz 

90* 

14.351 Hz 

14.247 Hz 

85.950 Hz 

85.934 Hz 
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From the above table it is observed that the hybrid model gives identical 
results for a rotation of 90‘ and no rotation of the material axes, 30* and 60* of the 
axes. The displacement method however gives results that vary, even though the 
moduli Ej and E 2 are equal. 

The first two mode shapes for the anisotropic cantilever, obtained using 
both the models do not vary much as shown in Fig. 15 (even for the case when the 
the difference between the eigen values is a maximum, viz., for 0 = 45°). 

Quadratic Element Results 

The natural frequencies and mode shapes of the isotropic and anisotropic 
cantilever beams are now calculated using an 8-noded finite element mesh with the 
number of elements varying from 3 to 20. The meshes used are the same as for the 
static case and are shown in Fig. 8. 

The normalized natural frequencies (coj anal /0)| f e and CO 2 anal f.e.) 
are plotted against the number of elements, and are shown in Fig. 16. Again, it is 
observed that the hybrid model converges to the analytical solution faster than the 
assumed-displacement method. The mode shapes however are very similar in both 
models, except for the maximum "amplitude" (when 10 quadratic elements are used) 
as shown in Fig. 17. 

The first two natural frequencies for various angles of rotation of the 
material axes, in an anisotropic cantilever beam, using both the displacement and 
hybrid approximations are tabulated in Table 4.7. The material properties and the 
material model assumed are the same as for the static anisotropic case, viz. 3 
independent constants in a crystal with cubic syngony, where 

Ej = E 2 = 1.9716 x 10 7 lb./in. 2 

v l2 = 0.2875 

G 12 = 5.4758 x 10 6 lb./in. 2 
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Table 4.7 


Natural frequencies for an anisotropic cantilever beam using 10 
quadratic elements for different material axes orientations. 


Orientation of the Axes 

0)j (hybrid) 

(Oj (displ.) 

©2 (hybrid) 

^(displ.) 

0* 

14.130 Hz 

14.197 Hz 

83.675 Hz 

84.682 Hz 

30* 

13.020 Hz 

14.061 Hz 

78.264 Hz 

83.995 Hz 

45* 

12.679 Hz 

13.058 Hz 

76.499 Hz 

78.911 Hz 

60* 

13.020 Hz 

13.707 Hz 

78.264 Hz 

82.211 Hz 

90* 

14.130 Hz 

13.276 Hz 

83.675 Hz 

80.023 Hz 


The results of the hybrid model are as expected, with frequencies falling 
as the angle of rotation is increased, reaching a minimum at a rotation of 45*, and then 
increasing symmetrically (since Ej = E 2 ) up to a rotation of 90*. 

The results of the displacement model do not show any symmetry about 
0 = 45* , and are not as invariant under a rotation of the axes. 

Again, there is not much difference in the mode shapes obtained from the 
hybrid and displacement models for the anisotropic cantilever, for 0 = 45°, as shown 
in Fig. 18. 

A Specific Numerical Example 

A tapered cantilever beam consisting of 3 crystals of the same material 
but with different orientations of the material axes is considered next, and analyzed 
for its displacements, stresses, natural frequencies and mode shapes. The beam is 
shown in Fig. 19. 

Since the Ni alloy for which experimental data was provided exhibits 
cubic syngony in its crystals, the same material properties as used in the previous 
sections are considered. Also, since of all the elements tested, the 8 noded 
hybrid-stress element gave the best results, it is used in the mesh shown in Fig. 20. 

The results are tabulated as follows: 
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Table 4.8 

Comparison of hybrid and displacement model results for a tapered, 
3-crystal anisotropic cantilever beam: 

Parameter Hybrid Method Displacement Method 


u tip 

a bending (max.) 
<c max 

CDi 

“>2 


0.0539 

1.632 x 10 5 psi 
1.643 x 10 5 psi 
30.198 Hz 
131.482 Hz 


0.0502" 

1.519 x 10 4 psi 
1.511 x 10 4 psi 
30.869 Hz 
143.976 Hz 


The displacements and stresses differ at most by about 6.5 percent, and 
the natural frequencies by even less. However, if the 30’ rotation is changed to a 60’ 
rotation and no rotation is changed to a 90’ rotation , all errors increase rapidly to 
almost 12-percent at most. 

The location of the points of maximum bending and shear stress are 
predicted accurately by both models, although the magnitudes predicted differ. The 
maximum bending stress is observed in elements 1,6 while the maximum shear stress 
is observed in elements 5,10 as the material axes in these two elements is rotated by 
45’ relative to the global frame. 

The mode shapes of the first two modes are shown in Fig. 21 but do not 
differ much, as in the previous cases discussed. 
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SUMMARY AND CONCLUSIONS 


In this work, a hybrid-stress finite element method is developed for 
equilibrium and vibration analysis of problems in two-dimensional anisotropic 
elasticity. 

A number of sample problems are solved using the hybrid-stress method 
and the standard displacement method and the results are compared. Emphasis is 
placed on a cantilever beam loaded in end shear because of its similarity to a turbine 
blade. 

It is observed that even for the isotropic case, the hybrid-stress model 
gives more accurate displacements, stresses and natural frequencies as compared to 
the results of the displacement method, although the variation between the two is not 
large. 

For anisotropic materials, especially when the material axes are rotated 
relative to the global axes, the hybrid-stress model is observed to be stable and 
invariant, while the displacement model shows some variation for pairs of rotations 
that are complementary when the two Young's moduli Ej and E 2 are equal. 

In the absence of analytical solutions for anisotropic cantilever beams, 
comparisons are made by increasing the number of elements and checking for 
convergence. The hybrid model behaves well even if the number of elements used is 
not too large, although if the degree of anisotropy is very high, e.g. Ej/E 2 = 10 4 , 
both models do not seem to converge rapidly. 

Work is now in progress to extend the finite element code to include three 

dimensional problems. The stress shape functions for 8-noded linear bricks and 
20-noded quadratic brick elements as proposed by Rubinstein, Punch and Atluri [6] 
will be implemented in the hybrid-stress finite element code and compared with the 
8-noded and 20-noded brick elements using a displacement approximation. 

Instead of using group theoretical methods to determine stable, invariant 
stress fields, complete stress polynomials may be chosen and the number of stress 
parameters reduced by forcing equilibrium and compatibility conditions to be 
satisfied. Although the algebra involved is tedious and the matrices are slightly 
stiffer, the element matrices will be very stable under rotation and the results thus 
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obtained may be compared with those obtained using the group theoretical stress 
polynomials. 

Once a three-dimensional finite element mesh has been generated for the 
anisotropic crystalline turbine blade, the above scheme may be used to analyze it 
statically and dynamically for stresses, displacements, natural frequencies and mode 
shapes. 

Another suggested development would be to investigate the use of 
triangular elements for two-dimensional hybrid-stress finite element analysis, and 
tetrahedrons for three-dimensional problems. 


APPENDIX A 


CALCULATION OF THE POLYNOMIAL STRESS FUNCTIONS 
FOR LINEAR (7 B) AND QUADRATIC (15 B) QUADRILATERAL 

ELEMENTS 


(a) linear quadrilateral elements: 

The stress polynomial is complete in linear terms and is given by 

o x = Bj + B 2 x + B 3 y 

a y = B 4 + B 5 x + B 6 y ( 1 ) 

X xy ~ + ^8 X + B 9 y 

In plane elasticity, the equilibrium conditions reduce to just two 
equations, viz. 


and 


dx 


Bt X y 

dy 


= 0 




(2) 


in the absence of body forces. 

Substituting (2) into (1) to ensure that the stress polynomials are 
equilibrated, we get 


Bg Hh 62 — 0 

% + % = 0 


(3) 


Eliminating B 8 and B 9 and renumbering the B's, we get the following 
equilibrated stress field for 4-noded quadrilateral elements; 
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( 4 ) 


a x = Bi + B 2 y + B 6 x 
a y = B 3 + B 4 x + Bjy 

X xy = ®7 X " 

b) quadratic quadrilateral elements: 

The stress polynomial is complete in cubic terms and is given by 

a x = Bj + B 2 x + B 3 y + 2 B 4 xy + B 5 x 2 + B 6 y 2 + B 7 x 3 + B 8 y 3 
+ 3B 9 xy 2 + 3B 1 Q x 2 y 

Oy = Bjj + B 12 X + B 13 y + 2B 14 xy + B 15 x 2 + B 16 y 2 + B 1 ? x 3 
+ B 18 y 3 + 3B 19 xy 2 + 36^^ 

T xy = B 21 + B 22 x + B 23 y + 2B 24 xy + B 25 x 2 + B 26 y 2 + B 27 x 3 

+ B 28 y 3 + 3B 29 xy 2 + 3B 3 Q x 2 y . (5) 

Substituting this stress field into the equilibrium equations ( 2 ), we get 

[B 2 + 2B 4 y + 2B 5 x + 3B 7 x 2 +3B 9 y 2 + 6 B 10 xy] + [B 23 + 

2B 24 x + 2B 26 y + 3B 28 y 2 + 6 B 29 xy + 3B 30 x 2 ] = 0 
and 

[B 13 + 2B 14 x + 2B 16 y + 3Bj 8 y 2 + 6 B 19 xy + 3B 20 x 2 ] + [B 22 + 
2B 24 y + 2 B 25 x + 3B 27 x 2 + 6 B 3 Qxy + 3B 29 y 2 ] = 0 

Equating coefficients of the polynomial terms separately to zero, we get 
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b 2 + b 23 — o : 

+ ®24 ” 0 ’ 

®4 + fi 26 = 0 ’ 


B 9 + ® * 

Bio + ®29 = ® * 

+ fi 30 = ® ’ 


®22 + ®13 = 0 ’» 

%4 + ^16 = ® > 

^14 + ^25 = o ; 


B 27 + B 20 = 0 ; 

29 + ®18 = ® * 

%0 + fii9 = 0 • 

(7) 


Eliminating B 22 through B 30 , B 16 , B 18 and B 19 and renumbering, we 
get the following 18B equilibrated stress field: 

o x = Bj + B 6 x + B 2 y + B g y 2 + 2B 9 xy + B^x 2 + B 13 x 3 + 

+ 3B 14 x 2 y + 3B 15 xy 2 + B 17 y 3 

c y = B 3 + B 4 x + B 7 y+ 2B n xy + B 12 x 2 + B l0 y 2 + 3B 13 xy 2 
+ B 14 y 3 + 3B 16 x 2 y 3 + B lg x 3 
T X y = B 5 - B 7 x - B 6 y - B^x 2 - B 9 y 2 - 2B 10 xy - SB^x^ 

- 3B 14 xy 2 - B 15 y 3 - B 16 x 3 (8) 

To reduce the number of B's still further, the stresses are allowed to 
satisfy the compatibility conditions necessary for the existence of a displacement 
field. 

In plane strain, there is just one compatibility relation, given by 

3% . 3% _d^xy 

3y2 3x 2 - 3x3y 

as expressed in terms of the strains e , e and y . Plane stress has additional 

x y xy 
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compatibility conditions, but these are in a sense more "restrictive" and not applied. 

For the most generalized anisotropic case, the stress-strain relations for 
plane elasticity are given by 


^x. ~ I a x ' v yx °y + ^xy.y ^yl 


ey = J_ [-v xy o x + o y + ri xy y Tx y ] 

Eyy 

1 r 

Y X y - t^x.xy a x + ^.xy CT y + X xyl 
^xy 


( 10 ) 


or, defining E = C a 

where C is the compliance matrix whose coefficients are 


1 1 1 

a ll “ p ’ a 22 ~ p ’ a 66 ~ q 

^xx ^yy ' J xy 


v yx ' v xy ^xy.x ^x.xy 

a 12 — g £ ’ a 16 £ £ 

c xx ^yy ^xx ^xy 


and 


a 26 = 


Oxy,y “Oy.xy 


yy 

^xy 


a n 

a 12 

a 16 


a 22 

a 26 » 

sym. 


a 66 


(ID 


T| x xy and T| yxy are coefficients of mutual influence of the second kind, and 
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Tlxy,x » Tlxy.y 816 coefficients of mutual influence of the first kind. 

To get the compatibility condition in terms of the stresses, we substitute 
equation (10) into (9) so that 

9 2 O x (PCy <Pl X y d 2 O x d 2 Oy 

+ 812 # + 816 + 812 U + 822 U 


9 2 x xy (Po x d 2 Oy ^ 2<c xy 

a 26 -j x 2 dxdy + fady + dxdy 


For an equilibrated stress field, 

d 2 o x d 2 Oy 23 2 x X y 
+ ^2~ _ dxdy 


(13) 


so that equation (12) becomes 


9 2 x 


xy 


d 2 a 


y , a 66 ^ 2 °y a 66 v v x 

+ , “*?- + < * , * + -’-V + < * b + -*->-5?- 


9 2 o v 


11 “ay 7 


+ a 16 


9 2 x 


xy 


d 2 X. 


d y 2 


^26 


xy 


dx 2 


= a 


16 


a 2g x 

dxdy H 


*26 


d 2 a y 
dxdy 


(14) 


Substituting the polynomial approximations for the equilibrated stresses 
from equation (8) , 
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Hjj (2Bg + 66 | 5 x + 6B 17 y) + 3j2 ( 2 Bj 2 ^ ®jg^) 

+ 2 (a 12 + ~ 2 ~ ) [26io + ^ 13 X + 6B 14 y] + a 16 (-2B 9 - 6B 14 x - 6B 15 y) 

+ a 26 (-2B 11 - 6B 13 y - 6B 16 x)= a 16 (2B 9 + 6B 14 x + 6B 15 y) 

+ ^26^^11 6Bi 3 y + 6B 16 x) (15) 

Equating the appropriate coefficients on both sides, we get 

a 11^8 + *22 ^12 + ( 2a 12 + a 66^ ®10 = 2a 16**9 + 2a 26®ll 

a ll Bj 5 + a^Bjg + ( 2a 12 + *66^ ®13 = 2a 16 ®14 + 2a 26**16 0^) 

a ll B17 + &22 ^16 + ( 2a 12 + a 66^14 = 2a 16 **15 + 2a 26 **13 


Eliminating B 12 ,B 17 and B lg from the above equations, and renumbering, we get a 
15 B cubic approximation for the intra-element stresses: 

C x = **l + ®6 X + B 2 y + V 2+ 2B 9 xy + B 10 x 2 

o / 3 2&26 3n „ ro 2. ( 2a 12 + a 66) 3l 

+ B 12 (x 3 + - — y 3 ) + B 13 [Sx-^y - : y 3 ] 


a ll 


a ll 


+ Bi4(3 X y 2 + ^y 3 ) B 15 y 3 


l ll 


a 


ll 


ll 


Cy = B 3 + B 4 x + B 7 y - — BgX 2 + B 10 [y z - 
3 a 22 


( 2a 12 + * 66 > 

*22 


x2 ] 


+ B 12 [3xy 2 - 


(2 a i2 + % 6 ) 

*22 


x 3 ] 


‘11 


22 


2a 


B 14 x 3 + B 15 (3x z y + 


15 


26 

*22 


x 3 ) 
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2a 


+ B 13 (y 3 + 


16 


*22 


x 3 ) + B 9 x 2 + B n (2xy + — — x 2 ) 
*22 *22 


T Jy = B 5 - 6,x - B 6 y - B„x 2 - B 9 y 2 - 2B 10 xy - 3B 12 x 2 y 

-3BUXX 2 - B 14 y 3 -B,jX 3 (17) 

For an isotropic material, 

a l6 = a 26 = 0 , and 2a 12 + a^ = 2a n 

a ll = a 22 ’ so that the stress field loses all dependence on the compliance 
constants and becomes simply 

o x = Bj + B 6 x + B 2 y + B g y 2 + 2B 9 xy + B 10 x 2 + B 12 x 3 
+ B 13 (3x2y - 2y 3 ) + 3B 14 xy 2 - B 15 y 3 
<y y = B 3 + B 4 x + B 7 y + 2B n xy - B g x 2 + B 10 (y 2 - 2x 2 ) 

+ B 12 (3xy 2 - 2x 3 ) + B 13 y 3 + 3B 15 x^ - B 14 x 3 
x xy = B 5 - B^ - B 6 y - B n x 2 - B 9 y 2 - 2B 10 xy - 3B 12 x2y 

- 3B 13 xy 2 - B 14 y 3 - B 15 x 3 (!8) 
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11 


( c ) 30 Elements , 44 Nodes 



( e ) 80 Elements , 123 Nodes 


Fig. 5 Various Meshes used with Linear 4-noded Elements 
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Fig. 6 Convergence of Beam Tip Displacements using 4-noded linear elements 
for Hybrid and Displacement models - Isotropic case 



Fig. 7 Convergence of maximum bending stress in the beam using 4-noded 
linear elements for Hybrid and Displacement models - Isotropic case 
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18 


3 



( a ) 3 Elements , 18 Nodes 



( b ) 5 Elements , 28 Nodes 




( d ) 20 Elements , 85 Nodes 


Fig. 8 Various Meshes used with 8-noded Quadratic Elements 
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0 5 10 15 20 25 


No. of Elements 

Fig. 9 Convergence of Beam Tip Displacement using 8-noded quadratic 
elements for Hybrid and Displacement models - Isotropic case 



5 10 15 20 25 


No. of Elements 

Fig. 10 Convergence of maximun bending stress in the beam using 8-noded 
elements for Hybrid and Displacement models - Isotropic case 
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i.oi 
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0 . 0 ' 
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Orientation of the material axes (in degrees) 


Fig. 11 Comparison of Beam Tip displacements for various rotations of the 
material axes using linear and quadratic Hybrid and Displacement elements 
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Fig. 1 3 Convergence of the first two Eigen values using 4-noded linear elements 
for Hybrid and Displacements models - Isotropic case 
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( a ) Mode 1 of the beam 



( b ) Mode 2 of the beam 


Fig. 14 The first two modes of vibration of an isotropic Cantilever Beam, obtained 
using 80 linear elements of the stress-hybrid and displacement model 
(not to scale) 






(b) Mode 2 of the beam 


Fig. 15 The first two modes for an anisotropic Cantilever Beam, (material axes 
rotated by 45 degrees), obtained using 40 linear elements of the stress- 
hybrid and displacement model (not to scale) 
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No. of Elements 
( a ) 1st Natural Frequency 



No. of Elements 


( b ) 2nd Natural Frequency 


Fig. 16 Convergence of the first two Eigen values using 8-noded quadratic 
elements for Hybrid and Displacement models - Isotrpic case 
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( a ) Mode 1 of the beam 



Fig. 17 The first two modes of an isotropic Cantilever Beam, using 10 quadratic 
elements of the stress-hybrid and displacement model (not to scale) 






( a ) Mode 1 of the beam 



Fig. 1 8 The first two mode shapes of an anisotropic Cantilever Beam (material 

axes roated by 45 degrees), obtained using 10 quadratic elements of the 
stress-hybrid and displacement model (not to scale) 
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150# 



Fig. 1 9 Tapered Cantilever Beam made up of 3 crystals of the same anisotropic 


X 
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( a ) Mode 2 of the beam 


ig. 21 The first two modes of the tapered Cantilever Beam made up of 3 
anisotropic crystals, obtained by using 10 quadratic elements of 
the stress-hybrid and displacement model (not to scale) 
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